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The magnetic energy-level diagrams for models of the Mni2 and V15 molecule are calculated 
using the Lanczos method with full orthogonalization and a Chebyshev-polynomial-based projector 
method. The effect of the Dzyaloshinskii-Moriya interaction on the appearance of energy-level 
repulsions and its relevance to the observation of steps in the time-dependent magnetization data 
is studied. We assess the usefulness of simplified models for the description of the zero-temperature 
, magnetization dynamics. 
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GO . I. INTRODUCTION 

a . 

Magnetic molecules such as Mni2 or V15 have attracted a lot of interest recently because these nanomagnets can be 
used to stu dy e .g. quantum (de)coherence, relaxation and tunneling of the magnetization on a nanoscale PISTOL 
T3i IE IE IE tlO, E3, O E3, EE EE E3, EE IS3, ISH As a result of the very weak intramolecular interactions 
between these molecules, experiments directly probe the magnetization dynamics of the individual molecules. In 
particular the adiabatic change of the magnetization at low-temperature is governed by the discrete energy-level 
structure Hill Ellll. 

The magnetic properties of molecules such as Mni2 or V15 are often studied by considering a simplified model for 
the magnetic energy levels for a specific spin multiplet, e.g. S=10 for Mni2 or S=3/2 for V15. However for these and 
other, similar, magnetic molecules that consist of several magnetic moments (12 in the case of Mni2, 15 in the case of 
V15), the reduction of the many-body Hamiltonian to an effective Hamiltonian for a specific spin multiplet is, except 
for the diagonal terms, non-trivial. 
'. Magnetic anisotropy, a result of the geometrical arrangement of the magnetic ions within a molecule of low symmetry, 
mixes states of different total spin and enforces a treatment of the full Hilbert space of the system. The dominant 
contribution to the mag netic anisotropy due to spin-orbit interactions is given by the Dzyaloshinskii-Moriya interaction 
(DMI) [27L I2H 12^. I30ll3ll |3| . In principle this interaction can change energy- level crossings into energy-level 
repulsions. The presence of the latter is essential to explain the adiabatic c hang es of the magnetization at the resonant 
fields in terms of the Landau-Zener-Stiickelberg (LZS) transition [23l l24l I25ll2r3|. Thus a minimal magnetic model 
Hamiltonian should contain (strong) Heisenber g in teractions, anisotropic interactions and a coupling to the applied 

magnetic field [is m m m m s m EE Ellil • 

In this paper we calculate the magnetic energy-level diagrams for models of the Mni 2 and V15 molecule using 
exact diagonalization techniques. We study the effect of the DMI on the appearance of energy- level repulsions that 
determine the adiabatic changes of the magnetization observed experimentally. In contrast to earlier work [37]>E1> the 
approach adopted in the present paper does not rely on perturbation theory. Instead we perform an exact numerical 
• i-H , diagonalization of the full Hamiltonian. 

■ As the quantum spin dynamics of these magnetic molecules is determined by the (tiny) level repulsions, a detailed 
knowledge of the low-lying energy levels scheme is necessary. In order to bridge the energy scales involved (e.g. from 
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500K, a typical energy scale for the interaction between individual magnetic ions, to w 10~ 2 — 10 _9 K, a typical energy 
scale for energy-level splittings), a calculation of the energy levels of these many-spin Hamiltonians has to be very 
accurate. We have tested many different standard algorithms to compute the low-lying states. For systems that are 
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FIG. 1: Schematic diagram of the dominant magnetic (Heisenberg) interactions of the Mni2 molecule. 
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FIG. 2: Schematic diagram of the magnetic interactions of the simplified model Q of the Mni2 molecule. 



too large to be solved by full exact diagonalization, we find that the Lanczos method with full orthogonalization and 
a Chebyshev-polynomial-based projector method can solve these rather large and difficult eigenvalue problems with 
sufficient accuracy. 

The paper is organized as follows. In Sec.[n|wc introduce the model Hamiltonians for the Mni2 and V15 molecules. 
In Sec. IIIII we briefly discuss the numerical algorithms that we use to compute the energy levels. Our results for the 
energy level schemes for Mni 2 and V15 are presented in Sec. II VI In Sec.0we analyze a reduced, 3-spin model for V15 
and determine the conditions on the DMI energy-level under which repulsions appear. Numerical calculations for the 
full V15 model confirm that these conditions are also relevant for the presence of energy-level repulsions in the V15 
model. 
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II. MODELS 



A. Manganese complex: Mni2 



In Fig. n we reproduce the schematic diagram of the dominant magnetic (Heisenberg) interactions of the Mn 12 
molecule (MniatCHaCOO^HaO^Oia^CHaCOOH^O). The four inner Mn+ 4 ions have a spin S = 3/2, the 
other eight Mn +3 ions have spin S = 2. The number of different spin states of this system is 4 4 x 5 s = 10 s . If the 
total magnetization is a conserved quantity, it can be used to block-diagonalize the Hamiltonian, allowing the study 
of models of this size [2H, However, to study the adiabatic change of magnetization, we have to treat all the 
states, and the dimension of the matrix become prohibitively large. Thus we need to simplify the model in order to 
reduce the dimension. A drastic reduction of the number of spin states can be achieved by assuming that the strong 
antiferromagnetic Heisenberg interaction (Ji) between an inner ion and its outer neighbor allows the replacement of 
the magnetic moment of an inner ion by an effective S=l/2 moment. The schematic diagram of this simplified (but 
still complicated) model is shown in Fig. El The number of different spin states of this model is 2 4 x 5 4 = 10 4 . In this 
paper we study the latter model. 

The Hamiltonian for the magnetic interactions of the simplified Mn 12 model can be written as [34^ 



4 2 4 8 

h = S 2 - 1 ) - Jl E S2i -! ■ s ^ - K * E (^) 2 + E D * J • t^- 1 x s ^ - E h ■ Si > w 

i=1 i=1 (id) 1=1 

where even (odd) numbered Si are the spin operators for the outer (inner) S = 2 (S = 1/2) spins. The first two 
terms describe the isotropic Heisenberg exchange between the spins. The third term describes the single-ion easy-axis 
anisotropy of S = 2 spins. In this paper we do not consider higher-order correction terms that restore the SU(2) 
symmetry [2!J, |3(j, |3l|, |44j . The fourth term represents the antisymmetric DMI in Mni2- The vector D 4 '- 7 determines 
the DMI between the i-th S — 1/2 spin and the j-th S = 2 spin. The last term describes the interaction of the spins 
with the external field h. Note that the factor g/iB is absorbed in our definition of h. 

The first three terms in Hamiltonian Q conserve the z-component of the total spin M z = Yli=i The DMI 
on the other hand mixes states with different total spin and also states with the same total spin. Hence, the DMI 
can change level crossings into level repulsions. Therefore, the presence of the DMI may be sufficient to explain the 
experimentally observed adiabatic changes of the magnetization. 

The four- fold rotational-reflection symmetry (S4) of the Mni2 molecule imposes some relations between the DM- 
vectors. It follows that there are only three independent DM-parameters: D x = D x ' s , D y = D^' 8 , and D z = D^ 8 , as 
indicated in Fig. [21 The above model satisfactorily describes a rather wide range of experimental data, such as the 
splitting of the neutron scattering peaks, results of EPR measurements and the temperature dependence of magnetic 
susceptibility |34|. The parameters of this model have been estimated by comparing experimental and theoretical 
data. In this paper we will use the parameter set B from Ref. pHioJ : J = 23.8K, J' = 79.2K, K z = 5.72K, D x = 22K, 
D v = 0, and D z = 10K. 

Although the amount of available data is not sufficient to fix all these parameters accurately, we expect that the 
general trends in the energy-level diagram will not change drastically if these parameters change relatively little. 



B. Vanadium complex: V15 



In Fig. |3| we show the schematic diagram of the dominant magnetic (Heisenberg) interactions of the V15 molecule 
(Kg[V{^As6042(H20)]-8H20). The magnetic structure consists of two hexagons with six S=l/2 spins each, enclosing a 
triangle with three S=l/2 spins. All dominant Heisenberg interactions are antiferromagnetic. The number of different 
spin states of this model is 2 15 = 32768. The minimal Hamiltonian for the magnetic interactions that incorporates 
the effects on magnetic anisotropy can be written as [23, H3, IS HU 



« = - E ■ Si + E D • & • s . E h - s t ■ ( 2 ) 



The various Heisenberg interactions Jjj are shown in Fig. |3 For simplicity, we assume that D JJ = for sites 
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and j except for bonds for which the Heisenberg exchange constant is J (see Fig. EJ 37, 41] . Rotations about 27r/3 
and 47r/3 around the axis perpendicular to and passing through the center of the hexagons leave the V15 complex 
invariant. This enforces constraints on the values of D 1 '- 7 l4li E3. In Sec. II VI we present results for several different 
sets of estimates for the model parameters of the V15 model |4l ll8l I37I El) . 
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FIG. 3: Schematic diagram of the magnetic interactions in model @ of the V15 molecule. 



III. NUMERICAL METHOD 

A theoretical description of quantum dynamical phenomena in the Mni2 and V15 nanomagnets requires detailed 
knowledge of their energy-level schemes. Disregarding the fascinating physics of the nanomagnets, the calculation 
of the eigenvalues of their model Hamiltonians is a challenging problem in its own right. Firstly, the (adiabatic) 
quantum dynamics of these systems is mainly determined by the (tiny) level repulsions. Therefore the calculation 
of the energy levels of these many-spin Hamiltonians has to be very accurate in order to bridge the energy scales 
involved (e.g. from 500K to « 1CP 9 K). Secondly, the level repulsions originate from the DMI that mix states with 
different magnetization. In principle, this prevents the use of the magnetization as a vehicle to block-diagonalize the 
Hamiltonian and effectively reduce the size of the matrices that have to be diagonalized. If a level repulsion involves 
states of significantly different magnetization (e.g. M z = — 10 and M z — 10) a perturbative calculation of the level 
splitting would require going to rather high order (at least 20), a cumbersome procedure. Therefore it is of interest 
to explore alternative routes to direct but accurate, brute-force diagonalization of the full model Hamiltonian. 

As a non-trivial set of reference data, we used the eigenvalues obtained by full diagonalization (using standard 
LAPACK algorithms) of the 10000 x 10000 matrix representing model Q 40]. For one set of model parameters, such 
a calculation takes about 2 hours of CPU time on an Athlon 1.8 GHz/1. 5Gb system. Clearly this is too slow if we 
want to compute the energy-level diagram as a function of the magnetic field h. In particular if we want to estimate 
the structure of the level splittings at the resonant fields we need the eigenvalues for many values of h. Furthermore, 
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in the case of V15 this calculation would take about 30 times longer and require about 15 Gb of memory which, for 
present-day computers, is too much to be of practical use. 

We have tested different standard algorithms to compute the low-lying eigenvalues of large matrices. The standard 
Lanczos method (including its conjugate gradient version) as well as the power method [4^. |48| either converge 
too slowly, lack the accuracy to resolve the (nearly)-degenerate eigenvalues, and sometimes even completely fail to 
correctly reproduce the low-lying part of the spectrum. This is not a surprise: by construction these methods work 
well if the ground state is not degenerate and there is little guarantuee that they will work if there are (nearly)- 
degenerate eigenvalues |47l l48j . In particular, the Lanczos procedure suffers from numerical instabilities due to the 
loss of orthogonalization of the Lanczos vectors [I?! |48| . It seems that model Hamiltonians for the nanoscale magnets 
provide a class of (complex Hermitian) eigenvalue problems that are hard to solve. 

Extensive tests lead us to the conclusion that only the Lanczos method with full orthogonalization (LFO) [4?l |4S| 
and a Chebyshev-polynomial-based projector method (CP) (see Appendix) can solve these rather large and difficult 
eigenvalue problems with sufficient accuracy. The former is significantly faster than the latter but using both gives 
extra confidence in the results. 

In the Lanczos method with full orthogonalization, each time a new Lanczos vector is generated we explicitly 
orthogonalize (to working precision) this vector to all, not just to the two previous, Lanczos vectors [47l |48|. With 
some minor modifications to restart the procedure when the Lanczos iteration terminate prematurely, after n steps 
this procedure tranforms nxn matrix H into a tri-diagonal matrix that is comparable in accuracy to the one obtained 
through Householder tri-diagonalization but offers no advantages 48] . In our case we are only interested in a few 
low-lying eigenstates of H. Thus we can exploit the fact that projection onto the (numerically exact) subspace of 
dimension k (k <^n), built by the Lanczos vectors will yield increasingly accurate estimates of the smallest (largest) 
eigenvalues and corresponding eigenvectors as k increases. 

In practice, to compute the M lowest energy levels, the LFO procedure is carried out as follows. 

• Perform a Lanczos step according to the standard procedure 

• Use the modified Gramm-Schmidt procedure to orthogonalize the new Lanczos vector with respect to all previous 

ones mm 

• Compute the matrix elements of the tridiagonal matrix 

• At regular intervals, diagonalize the tridiagonal matrix, compute the approximate eigenvectors <fi, y,i = 
(ipi\H\(fii) and Af — (fi\{H — [ii) 2 \(pi) for i = 1,...,M, and check if all are smaller that a specified 
threshold. If so, terminate the procedure (the exact eigenvalue Ei closest to ^ satisfies \Xi — A < E- L < /a; + A.;). 
If not, continue generating new Lanczos vectors, etc. 

IV. RESULTS 

A. Manganese complex: Mni2 

In Table [I] we present the numerical data for h = 0T and h = 5T, also obtained by LFO. The results obtained by 
full exact diagonalization (LAPACK), LFO and CP are the same to working precision (about 13 digits). In Fig. 0] 
we show the results for the lowest 21 energy levels of the Mni2 model as a function of the applied magnetic field as 
obtained by LFO. 

Although the total magnetization is not a good quantum number, we can label the various eigenstates by their 
(calculated) magnetization. For large fields and/or energies, eigenstates with total spin 8, 9 and 10 appear, as shown 
in Tabled In Fig. 0] eigenstates with \M Z \ ss 10(9) (within an error of about 10%) are represented by solid (dashed) 
lines (eigenstates with \M Z \ 8 appear for h > 4 but have been omitted for clarity). 

The standard S = 10 single-spin model for Mni 2 

U = -D{S z ) 2 -hS z , (3) 

is often used as a starting point to interpret experimental results 0, 0, Ull H2I Il3l [37j • The energy levels of this 
model exhibit crossings at the resonant fields h = ±Dn for n = —10, . . . , 10, in agreement with our numerical results 
for the more microscopic model JQ). For the parameter set B, we find that D w 0.55K, in good agreement with 
experiments 0,0- 

The single-spin model © commutes with the magnetization S z and therefore it only displays level crossings, no level 
repulsions. Adding an anisotropy term of the form S+ + St. only leads to level repulsions when the magnetization 
changes by 4, which does not agree with the observation of adiabatic changes of the magnetization for all h = 
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FIG. 4: The lowest 21 energy levels of the Mni2 model Q as a function of the applied magnetic field h. Solid lines: eigenstates 
with \M Z \ » 10; dashed lines: eigenstates with \M Z \ w 9. 



TABLE I: The 21 lowest eigenvalues Ei and total spin Si of the corresponding eigenstates of the Mni2 model (0 for two 
values of the external applied field h along the z-axis. The distance between Ei and the exact eigenvalue closest to Ei is 
At = (tfi\(H - Eif\(pi) 1/2 < 10~ 10 for i = 1, ... ,7. These calculations took about 20 minutes on an Athlon 1.8 GHz/1. 5Gb 
system, using 1000 fully orthogonal Lanczos vectors. 



i 


Ei (h = 0) 


Si(h = 0) 


Ei (h = 5T) 


Si (A = . 





-815.1971469173 


9.91 


-881.7827744750 


9.92 


1 


-815.1971469173 


9.91 


-860.4928253394 


9.93 


2 


-800.5810020061 


9.91 


-840.8569089483 


9.92 


3 


-800.5810020061 


9.91 


-822.8556918884 


9.92 


4 


-787.6124037484 


9.91 


-815.4339009404 


8.93 


5 


-787.6124037482 


9.91 


-811.0766283789 


8.93 


6 


-776.2715579413 


9.90 


-806.4609011890 


9.90 


7 


-776.2715579281 


9.90 


-797.9409264313 


8.94 


8 


-766.5314713958 


9.90 


-794.1268159385 


8.93 


9 


-766.5314702412 


9.90 


-791.6387794071 


9.90 


10 


-758.3618785887 


9.89 


-781.8373616760 


8.93 


11 


-758.3618126323 


9.89 


-778.4824830935 


8.96 


12 


-755.6412882369 


8.92 


-778.3500886860 


9.85 


13 


-755.6412882368 


8.92 


-776.4751565103 


8.93 


14 


-751.7362729420 


9.88 


-767.0677893890 


8.93 


15 


-751.7337526641 


9.88 


-766.5785427469 


9.87 


16 


-751.2349837637 


8.91 


-764.0838038821 


8.92 


17 


-751.2349837632 


8.91 


-761.4314952668 


8.76 


18 


-746.6655233754 


9.87 


-756.2910279030 


9.87 


19 


-746.6082906321 


9.87 


-753.5740765004 


8.92 


20 


-744.8208087762 


8.92 


-752.7461619357 


8.08 



nD 0,0,0,^3- I n contrast, for the DMI the Hamiltonian has nonzero matrix elements for the pairs of states \S, S z ) 
and \S ± 1, S z ± 1), but zero matrix elements for levels with the same value of the total spin. 

In Fig. for some values of h, level repulsions are present. However, these are due to the fitting procedure 
used to plot the data and the number of /i-values used (100) and disappear by using a higher resolution in h-ficlds 
(results not shown). Thus these splittings have no physical meaning. For the Mn 12 system, the energy splittings 
at low field are extremely small. Their calculation requires extended-precision (128-bit) arithmetic jijj. Therefore, 
to study the structure of the energy-level diagram in more detail we concentrate on the transitions at h rs 3.4T 
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FIG. 5: Left: The lowest 8 energy levels of Vis model (|5J with model parameters taken from Ref. |37| (VsetA) as a function 
of the applied magnetic field h parallel to the z-axis. Right: Detailed view of the four lowest energy levels at h ~ 0. 



TABLE II: The eight lowest eigenvalues Ei and total spin Si of the corresponding eigenstates of the Vis model © with model 
parameters taken from Ref. 37] (VsetA) for two values of the external applied field h parallel the z-axis. The distance between 
Ei and the exact eigenvalue closest to Ei is Aj = {ipi\(H — Ei) 2 \ipi) x ^ 2 < xl0~ 9 for i = I, ... ,7. These calculations took less 
than 20 minutes on a Cray SV1 computer, using 521 fully orthogonal Lanczos vectors. 



i 


Ei(h = 0) 


Si(h = 0) 


Ei(h = 4T) 


Si(h = 4T) 





-3679.53623744 


0.51 


-3683.51181131 


1.50 


1 


-3679.53623744 


0.51 


-3682.21997451 


0.51 


2 


-3679.52777009 


0.51 


-3682.18488706 


0.53 


3 


-3679.52777009 


0.51 


-3678.11784886 


1.50 


4 


-3675.42943612 


1.50 


-3676.84225573 


0.52 


5 


-3675.42943612 


1.50 


-3676.83951808 


0.51 


6 


-3675.42325141 


1.50 


-3672.74011178 


1.50 


7 


-3675.42325141 


1.50 


-3667.37940477 


1.50 



[M z k, —10 — > M z w 4)and h rs 3.9T fM z : sa —10 — ► M z sa 3) for which adiabatic changes of the magnetization have 
been observed in experiments 0, 0, Ull E3 • From experiments one finds that the magnitude of these splittings is of 
the order of 10 nK 45]. Extensive calculations lead us to the conclusion that the energy splitting at these resonant 
fields is smaller than 10 _6 K. Adding an extra transverse field by tilting the h-field by 5 degrees does not change 
this conclusion. Thus, it is clear that within the (very high) resolution in the h-field and 13-digit precision of the 
calculation, there is no compelling evidence that the DMI gives rise to a level repulsion, at least not for the choice of 
model parameters (set B, see above) considered here. The algorithms developed for the work presented in this paper 
can be used for 33-digit calculations without modification and we leave the calculation of the splittings for future 
work. 



B. Vanadium complex: Vis 



For the model parameters given in Ref. [33, J = -800, Ji = J' = -54. 4K, and J 2 = J" = -160K, J 3 = 
J4 = J5 = Jq = and in the absence of the DMI, we find for the energy gap between the ground state and 
the first excited state at h = a value of 4.12478K, in perfect agreement with Ref. [37]]. Following Ref. [4^ we 
take for the DMI parameters D^ 2 = Dy' 2 = D\ 2 — 40K, which is approximately 5% of the largest Heisenberg 
coupling. Using the rotational symmetry of the hexagon we have D^ A — 14.641K, Dy 4 = — 54.641K, — 40K and 
L>5< 6 = -54.641K, D^ 6 = 14.641K, D 5 Z ' 6 = 40K. As the two hexagons are not equivalent we cannot use symmetry 
to reduce the number of free parameters. For simplicity, we assume that the (x,y) positions of the spins on the 



8 



-3666 
-3668 
-3670 
-3672 
g -3674 
w -3676 
-3678 
-3680 
-3682 
-3684 



0.5 1 



1.5 2 2.5 
h[T] 



3 3.5 4 



-3679.46 
-3679.48 

-3679.5 
-3679.52 
-3679.54 
-3679.56 
-3679.58 

-3679.6 
-3679.62 



0.02 0.04 0.06 0.08 
h[T] 



0.1 



FIG. 6: Left: The lowest 8 energy levels of V15 model J5J with model parameters taken from Ref. |41| (VsetB) as a function 
of the applied magnetic field h parallel to the z-axis. Right: Detailed view of the four lowest energy levels at h ~ 0. 



lower hexagons differ from those on the upper hexagon by a rotation about 7r/3. This yields for the remaining model 
parameters Df' 11 = -14.641K, Dy ' 11 = 54.641K, D 1 ®' 11 = 40K, D* 2 ' 13 = -40K, D^' 13 = -40K, D] 2 ^ = 40K, 
and D^ 15 = 54.641K, D* 4 * 15 = -14.641K, L>] 4 < 15 = 40K. We will refer to this choice as VsetA. In Fig. we show 
the results for the eight lowest energy levels of V15 model as a function of the applied magnetic field along the 
z-axis, using the parameters VsetA. 

From Table ITT1 we see that for zero held, the DMI splits the doubly-degenerate doublet of S — 1/2 states into 
two doublets of S = 1/2 states. The difference in energy between the doubly-degenerate, first excited states and 
the two-fold degenerate ground states is due to the DMI and, for the parameters VsetA, has a value of w 0.0085K, 
much smaller than the experimental estimate w 0.05K 22], but of the same order of magnitude as the values cited in 
Ref. 0]. The next four higher levels are S = 3/2 states. The energy-level split ting between the S — 3/2 and S =1/2 
states is « 4. IK, in reasonable agreement with the experimental value 3.7K [45|. 

Following Ref. 0, we take J = -800, J\ = J' = -225K, J 2 = J" = -350K, and J 3 = J 4 = J5 = J 6 = 0. In the 
absence of a DMI, we find that the energy gap between the four-fold degenerate ground state and the first excited state 
is 3.61K, in ful l ag reement with the result of Ref. ^lj. Note that this value of the gap is fairly close to the experimental 
value of 3.7K 0. Taking for the non-zero DMIs D l x ' 2 = D 1 ^ 15 = 25K, L> 3 > 4 = D 5S = D io,n = £,12,13 = _ 12 .5K, 
D 3 ' 4 = —Dy' 6 = —Dy ' 11 = D 4 2 ' 1 3 = — 21.5K, our calculation for the splitting between the two doubly-degenerate 
S=l/2 levels yields 0.0037K, about a factor of two larger than the value cited in Ref. 0]. F° r the energy splitting 
between the S — 1/2 and S — 3/2 levels we obtain 3.616K instead of the value 3.618K given in Ref. |41l. These 
differences seem to suggest that a perturbation approach for the DMI has to be applied with great care |4y| . In Fig. 
we show the results for J = —800, J\ = J' = — 225K, and J 2 = J" = — 350K |4J and the same DMI parameters as 
in VsetA (which we will refer to as VsetB). 

For the energy gap at zero field, we find 4. IK and 3.61K for VsetA and VsetB respectively whereas the experi- 
mental estimate is 3.7K ^5|. The transition between the states 1 1/2, 1/2) and |3/2,3/2) takes place at h w 2.8T and 
h w 3.0T respectively, also in good agreement with the experimental value 2.8T. 

The most advanced estimation of the model parameters VsetC is given in Ref. [lg- Taking J = —809, J' = — 120K, 
J" = 120K, J x = -30K, J 2 = -122K, J 3 = -3K, Ja = -UK, J 5 = -3K, J 6 = -2K (see Table I in Ref. [3) yields 
an energy gap of 4.915K, in agreement with Ref. 18]. At h w 3.6T, the S = 1/2 and S = 3/2 states mix, a level 
repulsion appears and the adiabatic change of the magnetization from M w 1/2 to M k 3/2 gives rise to a step in 
the magnetization versus (time-dependent) h-field. Although the qualitative features of the energy-level diagram for 
VsetC also agree with what one would expect on the basis of experiments, the field at which the states 1 1/2, 1/2) 
and 1 3/2, 3/2) cross, h « 3.6T, does not compare well to the experimental estimate h ~ 2.8T. 

On a coarse scale, the level diagrams for VsetA, VsetB and VsetC are all similar and also resemble those of 
Ref. However, on a finer /i-scale a new feature appears (see right panel of Figs. EI EI andEJ): the field at which 
the energy difference between the second and third level reaches a minimum is no longer at h = 0. In other words, 
in the presence of the DMI, the adiabatic transition between the states 1 1/2, —1/2) and 1 1/2 , 1/2) does not occur. As 
we show in the next section, this seems to be a generic feature of the DMI in models of V15. 
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FIG. 7: Left: The lowest 8 energy levels of V15 model © with model parameters taken from Ref. (VsetC) as a function 
of the applied magnetic field h parallel the z-axis. Right: Detailed view of the four lowest energy levels at h ~ 0. 



V. DISCUSSION 



As shown above, the effect of the DMI on the energy-level diagram is much larger for the V15 model than it is 
for the Mni2 model. Therefore, to study the possibility of using simplified models for capturing the essential time- 
dependent magnetization dynamics, we will focus on models for the V15 molecule which is somewhat easier to treat 
numerically. Qualitatively the energy-level scheme for the eight lowest energy levels of the V15 models considered in 
Sec. IIVI closely resembles the energy- level diagram of a reduced, anisotropic model of three 5 = 1/2 spins described 
by the Hamiltonian 0, El 

3 

H = —J (Si ■ S 2 + S 2 • S 3 + Si ■ S 3 ) + D 1 ' 2 ■ [Si x S 2 ] + D 2 ' 3 ■ [S 2 x S 3 ] + D 1 3 • [S x x S 3 ] - ^ h • (4) 

i=l 

In the absence of the DMI, fitting the energy-level diagram of model (0J to exprimental data yields J ss — 2.5K [l5| . 
We use this estimate to fix J in our numerical calculations. The number of free parameters can be reduced further 
by exploiting the rotational symmetry of the triangle. We have D x ' 2 — D x , D y ' 2 — D y , D 2 ' 3 = — (^/3D X + D y )/2, 
D^ 3 = -(D x + V3D y )/2, D\/ = (V3D X - D y )/2, and D\ 2 = D 2 / = 0_ 3 = D z . The numerical results presented 
in this paper have been obtained for D x = D y = D z = 0.1K. In Ref. [23 the DMI vector is taken parallel to the 
y-axis at all the bonds and the field is applied along to the z-axis. This case corresponds to the case with only D z in 
the present model with the field applied in the x-direction. In this case the gap opens symmetrically with field .22|. 
However, as we show in this paper, the structure of the gap depends on the direction of the field. 

In Fig. we present results for the eight lowest energy levels of the three-spin model Q as a function of the applied 
magnetic field along the z-axis. Qualitatively it agrees with the level diagram of the full Vis model with parameters 
VsetA. The effect of the DMI is two-fold: as expected it lifts degeneracies but it may also shift the position of the 
resonant points in a non-trivial manner. A similar effect was also found in the full model calculations (see Sec. II V|) . 

The butterfly hysteresis loop observed in time-resolved magnetization measurements has been interpreted in terms 
of combination of a LZS transition at zero field and spin-phonon coupling |l5|,|22j. Here it should be noted that unless 
the field is applied in a special direction (x or y direction in this case), the set of avoided level crossings is no longer 
symmetric with respect to the field. Indeed, a closer look at the level diagram (see left picture in Fig.[SJ) reveals that 
the mimimum energy difference between the two pairs of levels does not occur at zero field but at h « 0.05T. This 
implies that the LZS transition from |l/2, —1/2) to the 1 1/2, 1/2) level does not take place at h = but at h « 0.05T. 
The minimum energy splitting between the first and second level (counting states starting from the ground state) 
also depends on the direction of the field. For the model parameters used in our calculations, it increases from 0.05T 
for h parallel to the z-axis to 0.12T for h parallel to the x-axis (results not shown). The fact that the DMI not only 
lifts the denegeneracy but, depending on the direction of the field with respect to the symmetry axis, also shifts the 
resonant point away from h = seems to be a generic feature. 

Summarizing: Our numerical data for the parameters VsetA, VsetB, and VsetC suggest that the three-spin 
model reproduces the main features of the full V15 model. The presence of the DMI allows for adiabatic changes of 
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FIG. 8: Left: The eight lowest energy levels of V15 model @ asa function of the applied magnetic field h parallel to the z-axis. 
Right: Detailed view of the four lowest energy levels at h « 0. Note that the energy- level splitting between the second and 
third level reaches a minimum at h « 0.05T, not at h = 0. 



the magnetization but, according to our calculations, the value of the resonant field for the |l/2, —1/2) to 1 1/2 , 1/2) 
transition changes with the direction of the magnetic field. This change (by a factor of two at least) should lead to 
observable changes in the hysteresis loops but has not been seen in experiment 45]. Therefore, although the DMI 
causes the avoided level crossing structure, it is anisotropic with respect to the direction of the field. Within the three 
spin model we have studied the effects of higher-order correction terms that restore the SU(2) symmetry |2^ . l8fll8lll43 |. 
and found that it has no essential effect on the low energy degenerate doublets while it causes the four S = 3/2 levels 
to be degenerate at h = 0. In experiments only weak directional dependence was found. Thus, another type of 
mechanism for the gap such as hyper- fine interaction, etc., is necessary and will be studied in the future. 
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Appendix: Projection method 

As an alternative to the Lanczos method with full orthogonalization, we have used a power method 0, ^3 based 
on the matrix exponential e~ tH [49j. Writing the random vector \I/(0) in terms of the (unknown) eigenvectors {<f>i} of 
H, we find 



V(t) = e- tE ° [0 o (0o|*(O))+e-*( Bl ^V 1 (0 1 |vl/(O))+ e - t ^- Bc ')0 2 (0 2 |v|/(O)) + ...] , (5) 

showing lim^oo $(t)/||$(t)|| oc <fi if (0q|\I>(O)) 7^ 0. In this naive matrix-exponential version of the power method, 
convergence to the lowest eigenstate is exponential in t if E\ > Eq . 

The case of degenerate (E = E\ = ...) or very close (E m E\ m ...) eigenvalues can be solved rather easily by 
applying the projector to a subspace instead of a single vector, in combination with diagonalization of e tH within this 
subspace First we fix the dimension k of the subspace by taking k equal or larger than the desired number of 
distinct eigenvalues. The projection parameter t should be as large as possible but nevertheless sufficiently small so 
that at least the first k terms survive one projection step. Then we generate a set of random initial vectors ^(O) for 
i = 1, . . . , k and set the projection count n to zero. We compute the k lowest eigenstates by the following algorithm [4fll | 

• Perform a projection step ^>i((n + l)t) = e" tH ^ i{nt) for i = 1, . . . , k. 

• Compute the k x k matrices. A = (*. 4 ((n + l)t)\e tH |*j((n + l)t)) = (*,-((n + l)t)\^i(nt)) and B = (*i((n + 
l)t)|4'i(n + l)i)). Note that A is hermitian and B is positive definite. 
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• Determine the unitary transformation U that solves the k x k generalized eigenvalue problem Ax — XBx. Recall 
that k is small. 

• Compute %({n + l)t) = £* =1 UijVj((n + l)t) for i = 1, . . . , k. 

• Set *i((n + l)t) = #<((n + l)t) for i = 1, . . . , k. 

• Compute ^ = l)t)|fl"|*i((n + 1)*)) and check if Af = l)t)|(il - /ij) 2 |*i((n + 1)*)) is smaller 
than a specified threshold for i = 1, . . . , &. If yes, terminate the calculation. If no, increase n by one and repeat 
the procedure. 

We calculate e~ tH & by using the Chebyshev polynomial expansion method [HoL IBH l52l IH3. l54| . First we compute 
an upperbound R of the spectral radius of H (i.e., \\H\\ < R) by repeatedly using the triangle inequality 0. From 
this point on we use the "normalized" matrix H = (2H/R — l)/2. The eigenvalues of the hermitian matrix H are 
real and lie in the interval [—1, 1] |47l |48|. Expanding the initial value ^(0) in the (unknown) eigenvectors 4>j of H 
(or H) we find 

*(t) = e -* H *( ) = e ***(0) = 5>'^<&|*(0)), (6) 

3 

where z = —tR. We find the Chebyshev polynomial expansion of ^(t) by computing the Fourier coefficients of the 

function e 2COs9 [55|. Alternatively, since — 1 < Ej < 1, we can use the expansion e zEj = 4(z)+2Xim=i^ni( z )^»(^) 
where I m {z) is the modified Bessel function of integer order m jH^ to write Eq. © as 



*(*) = 



I {z)I + 2j2 I m(z)T m (H) 



*(0). (7) 



Here, / is the identity matrix and T m (H) is the matrix- valued Chebyshev polynomial defined by the recursion relations 

T„(fl)¥(0) = #(0), ri(#)*(0)=fl-*(0), (8) 

and 

T m+ i(H)*(0) = 2HT m (H)V(0) - T m _i(fl)*(0) , (9) 

for m > 1. In practice we will sum only contributions with m < M where M is choosen such that for all 
m > M, \I m (z) / Iq(z)\ is zero to machine precision. Then it is not difficult to show that \\e~ tH / Io(z) — I — 
2^ m=1 [I m {z) / Io{z)]T m (H)\\ is zero to machine precision too (instead of e~ tH we can equally well use e~ tH /I (z) as 
the projector). 

Using the downward recursion relation of the modified Bessel functions, we can compute K Bessel functions to 
machine precision using only of the order of K arithmetic operations [5^, Ei|. A calculation of the first 20000 
modified Bessel functions takes less than 1 second on a Pentium III 600 MHz mobile processor, using 14-15 digit 
arithmetic. Hence this part of a calculation is a negligible fraction of the total computational work for solving the 
eigenvalue problem. Performing one projection step with e~ tH amounts to repeatedly using recursion to obtain 
T m (B)^ (0) for k = 2, . . . , M, multiply the elements of this vector by I m (z) and add all contributions. 
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